Multifractal analysis of normal RR heart-interbeat signals in power spectra ranges. 
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Power spectral density is an accepted measure of heart rate variability Two estimators of mul- 
tifractal properties: Wavelet Transform Modulus Maxima and Multifractal Detrended Fluctua- 
tion Analysis are used to investigate multifractal properties for the three strongly physiologically 
grounded components of power spectra: low frequency (LF), very low frequency (VLF) and ultra 
low frequency (ULV). Circadian rhythm changes are examined by discrimination of daily activity 
from nocturnal rest. Investigations consider normal sinus rhythms of healthy 39 subjects which are 
grouped in two sets: 5-hour wake series and 5-hour sleep series. Qualitative arguments are provided 
to conjecture the presence of stochastic persistence in LF range, loss of heart rate variability during 
night in VLF range and its increase in ULF. 

PACS numbers: 87.19.Hh, 87.80.Vt, 89.75.Fb, 05.45.Pq 
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I. INTRODUCTION 

Power spectral analysis of normal RR- series of heart 
rate has gained popularity as a simple, non-invasive mea- 
surement capable of assessing dynamic changes in the au- 
tonomic nervous system's control of heart rate [7]. Power 
spectrum means frequency analysis to identify superim- 
posed oscillations which contribute to variations in heart 
rate. It has been observed [5] that in the very low fre- 
quency range, from 0.0003 to 0.04 Hz, the power spec- 
tral density decays as 1//^. This property provides evi- 
dence of long range dependences in the data which, when 
viewed within the right range of scales, leads to the no- 
tion of fractality or self-similarity: statistical properties 
measured at some time scale look just like an appropri- 
ately scaled same statistical properties measured over a 
different time scale. 

However long range dependences leave the small-time 
scale behavior (high frequency and low frequency compo- 
nents of the spectrum) essentially unspecified. In order 
to be capable of accounting for small time scales we need 
methods which deal with regularity, it is with differentia- 
bility, of stochastic processes. Traditionally to quantify 
and qualify local irregularities the spectrum of Holder ex- 
ponents, called multifractal spectrum is indicated @, 3- 
The practical estimations of exponents is performed in 
three steps. First, from the observed time series X(t) 
one computes multiresolution coefficients Tx(a,t), quan- 
tities that depent jointly on time t and an analysis scale 
a. Second, one constructs the partition function P(g, a) 
from these coefficients weighed by q powers. Third, one 
measures the slope T{q) in log a of versus logP(g', a) dia- 
gram. The multifractality of X{t) means that r(g) ^ qH. 
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The aims of the following are: 

— To provide links between properties learned from 
multifractal studies and characterization of the heart con- 
trol system obtained by power spectral density analysis 
in physiologically grounded frequency intervals. Espe- 
cially, daily activity will be compared to nocturnal rest. 
By separating wake from sleep human stages we want 
to minimize effects related to nonstationarity in a signal. 
So far the problem whether the multifractal properties 
depend on a man activity is discussed and contradicting 
conjectures are suggested [E", 'g^. 

— To present in a compact and uniform way the nu- 
merical techniques used in discovering multifactal prop- 
erties. We discuss results arisen from two most popular 
methods: Wavelet Transform Modulus Maxima [7] and 
Multifractal Detrended Fluctuation Analysis [8]. More- 
over, following [9], we consider the linearization effect 
in scaling of a structure function r{q) by searching for 
points where a multifractal spectrum concentrates. 

The estimates are applied to two types of series: to a 
normal RR-signal extracted from ECG recordings and to 
the series obtained by integration of a normal RR-signal. 
The multifractal properties are statistical features. They 
are obtained numerically following heuristic procedures. 
So that results could depend on details of both - prop- 
erties of series and numerical procedures. Therefore, for 
reliability of the discussed multifractality it should be 
important to confirm results in various ways [lo] . 

The paper is organized as follows. In Sec. II we give 
a brief physiological background to reasons of heart rate 
variability. In Sec. HI the multifractal formalism is intro- 
duced and discussed. Then, Sec. IV, we shortly describe 
data which is used in our investigations. Sections V and 
VI contain the results: the group partition functions to- 
gether with variety of linear fits which can be found to 
these functions (Sec. V) and the multifractal spectra es- 
timated in physiologically grounded time intervals. The 
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FIG. 1: EGG diagram for a single heart period. Subsequent 
elements corresponds to atrial depolarization: P-wave, ven- 
tricular depolarization: QRS-complex and ventricular repo- 
larization: T-wave. 




FIG. 2: A series of RR normal heart beat intervals (green) 
and the same data shuffled at random. The series is 5-hour 
record of daily activity of a healthy subject (nsr_gda03_wake 

M) 



last section, Sec. VII, presents suggestions for further 
work. We believe that our considerations, though con- 
centrated mostly on the purely descriptive aspect of mul- 
tiscaling, are important in putting forward the research 
aimed on reasons for complex dynamics that is observed 
by the heart rate. 



II. HEART RATE VARIABILITY BY POWER 
SPECTRUM ANALYSIS 

The timing of a normal heart beat consists of three 
elements easily recognizable in an electrocardiogram 
(EGG), see Fig. [H In a normal heart beat the following 
sequence of events must occur: P-wave (describing depo- 
larization of atria), QRS-complex (associated to contrac- 
tion of ventricles) and then T-wave (showing relaxation 
of ventricles after the contraction). The so-called nor- 
mal sinus rhythm means that each P-wave is followed by 
QRS-complex and the rate is characterized by 60 to 100 
beats per minute with 10% variation [1]. The time in- 
tervals between successive R-peaks of QRS-complex are 
called RR-series. Heart rate variability refers to increases 
and decreases over time in the RR interval. In Fig. [2] a 
RR-series is compared to the same data but shuffled at 
random. An-eye inspection hints that there must be a 
relation between consecutive heart beats. 

A heart beat is generated by the autonomous pace- 
maker, called sino- atrial node, located at the top of the 



right atrium. However, the frequency is mediated by 
the neural (autonomic nervous system — ANS) and hor- 
monal Jrennin - angiotensin - aldosterone system) activ- 
ity [l|, [ll] . This activity is controlled by large number of 
different feedbacks acting simultaneously. For example, 
both the blood pressure as well as changes in the blood 
pressure stimulates baroreceptors [l2| , and baroreceptors 
generate a signal to ANS. 

ANS consists of two parts which act antagonistically — 
the sympathetic subsystem activity accelerates the heart 
rate and enhances the heart muscle contractility while the 
parasympathetic subsystem activity elongates the inter- 
beat intervals. This permanent conflict is supposed to 
be the reason for complex fluctuations which the normal 
heart rate dis play s. It is a widely accepted conjecture 
[l|, [ll L [12I, [13I, | l4 [ [15I, fl6ll , strongly supported by physi- 
cists [gI, [13, [m, [isl, [20, UM ^^^^ variability of the heart rate 
reveals the delicate balance between the two branches of 
ANS and even more — 1// scaling in healthy heart rate 
requires the existence of this balance [22| . 

A recent introductory to the methodology of the heart 
rate variability can be found in [15] and discussion of 
autonomic system influence in [3]. Here we concen- 
trate on the so-called frequency-domain measures re- 
sulted from analysis of power spectral density of a RR 
tachogram [1]. To avoid nonstationarities short-time 
recordings (several minutes) are used. Three main spec- 
tral components are distinguished in this spectrum: very 
low frequency: 0.003 < / < 0.04 Hz (VLF), low fre- 
quency 0.04 < / < 0.15 Hz (LF) and high frequency: 
0.15 < / < 0.4 Hz (HE) components. Longer record- 
ing (hours) spectral analysis provides an additional range 
/ < 0.003 Hz termed ultra low frequency component 
(ULE). 

The physiological explanation of ULE and VLF is less 
defined than other components. The power spectral den- 
sity for these components is of power-low form 1//^. The 
value of /3 can discriminate healthy heart rate from a 
rate with heart failure [2]. Moreover, reduced levels of 
the spectral power have been identified as predictions 
of all-cause mortality (see [l6| for discussion). There is 
some evidence that it is influenced by thermoregulatory 
mechanisms and humoral rennin-angiotensin system, but 
physical activity also notably influences to the power [26]. 

The higher frequencies — HE component, are mod- 
ulated predominantly by parasympathetic activity. A 
sino-atrial node responds to this activity at the respi- 
ratory frequency. LF component is modulated by both 
sympathetic and parasympathetic neural influence. The 
ratio of the spectral power LF to HE is used as a mea- 
sure of sympathovagal balance. However this measure is 
found to be controversial [l5^. 

It is known that multifractal analysis can serve indica- 
tors of a state of ANS. The definite reduction of multi- 
fractal properties, namely the transition from multifrac- 
tal to monofractal, is observed in multifractal picture of 
EGG of patients with congestive heart failure [13, fisj . 
The modifications in multifractal spectra are observed 
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in case of other heart diseases [TqI, [20|, [23|. Fohowing 
those investigations some heart rate variabihty measures 
have been proposed, e.g , spectrum width, spectrum lo- 
cation, curvature of a spectrum. Especiahy, it has been 
observed that a multifractal spectrum is distinct from a 
monofractal plot by existence of points (more than one) 
at which the spectrum is concentrated [20] . Speculations 
are continued that multifractal analysis offers potentially 
promising tools for heart rate variability assessments, but 
standards are lacking. 

III. MULTIFRACTAL APPROACH 

Multifractal analysis of time series is strongly math- 
ematically grounded, see for example 0, 1^. Roughly 
speaking by multifractal study one wants to express rela- 
tions between points X{t') and a given point where 
t' e J^{t) belongs to some neighborhood of t, by compar- 
ing a part of series to a polynomial with a singularity 
exponent h{t): 

\X{t')-X{t)\^ K\t' -t\^^'^ for t'eM{t) (1) 

The singularity exponent h is often called Holder expo- 
nent. Typically a process X{t) will posses many different 
singularities when moving t. The multifractal spectrum 
measures the frequency of occurrence of a given singular- 
ity exponent. Since singularity values change along the 
series vividly then the probability of occurrence is mea- 
sured by the Hausdorff dimension of the subset of time 
where the same singularity exponent value is accounted: 

D{h) = dimnit : h{t) = h} (2) 

However transferring these ideas to numerical estimates 
of a series which arises from some natural phenomena 
is not obvious. Both tasks: extracting singularity expo- 
nents and evaluating the Hausdorff dimension, are nu- 
merically difficult and calculations cannot be performed 
automatically. 

When one deals with discrete time series then special 
difficulty relates with determination of the neighborhood 
JV{t) . In case of continuous data the fractal sin- 
gularity is estimated with respects to the small enough 
\t — t'\ interval. But how to transfer this condition to a 
discrete signal? At the smallest possible interval, namely 
a unit interval, each two points can be connected by any 
function. If a regular function is chosen then the fractal 
dimension will be zero. Does it mean that the described 
procedure is useless? In Figl3]we present a piece of series 
consisting of a hundred points which are divided into two 
equal size boxes. The question is whether the neighbor- 
hood consisting of 50 points or more can be considered 
as small enough to provide estimates for roughness of a 
series. In our opinion any scale is appropriate to discuss 
the fractal properties of a discrete time series. Following 
this belief we will continue to call our search as calcula- 
tions of multifractal properties. 
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t+DO 1^0 150Q 

nuni-ber of a l>eart beat 

FIG. 3: WTMM vs MDFA methods in a simplest form. A 
signal is divided into boxes of size 50 heart beats. In each box, 
the oscillations of a signal are represented by two estimates: 
a wavelet coefficient (here, the Haar wavelet is plotted by a 
solid line) and a fluctuation from a main trend (a dashed line 
plots the linear trend for a given box). 



Both estimates: singularity exponent h and its fractal 
dimension D{h) can be found thanks to scaling properties 
of the partition function and their relation to a global 
characterization of a fractal object [3|, [2^. 

Basically, there are two approaches to construct the 
multiresolution coefficients: 

(1) aggregation of a signal within a i^-th box of a given 
size n, i.e., X(n, K) = ^ Z^ILi ^{^^ + '^)^ 

and 

(2) increment of a signal within a i^-th box of a given 
size n, i.e., X(n, k) = X{Kn + n) — X{Kn). 

By the wavelet technique and detrended fluctuation 
method one deals with both: increments between care- 
fully chosen averages are estimated. Having a signal 
divided into boxes of a given scale, see Fig. [3l the 
wavelet transformation in a i^-th box: T^{n^K) mea- 
sures the signal local oscillations inside a box. One can 
also estimate signal oscillations by measuring departures 
FpK (n, K) of the signal from a polynomial trend of 
order m found for points of i^-th box. Finally, the par- 
tition function for a given scale is defined as the sum of 
either selected — local absolute value maxima in case 
of wavelets, or all — in case of detrended fluctuations, 
oscillations weighted by q exponent. Both methods are 
heuristic and therefore below we give the precise mathe- 
matical formulas of estimates. 

Let us assume that a series {X{i)} is divided into boxes 
consisting of n points. Then the following calculations 
lead to the partition functions Z{n^q) and F{n^q), re- 
spectively. 

— In Wavelet Transform Modulus Maxima (WTMM) 
method a multiresolution coefficient is obtained as fol- 
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lows, 0: 



(3) 



where '^^(n, i) is the mother wavelet transformed to scale 
n and shifted to i^th box. Then the partition function 
q) is: 



Z{n,q) 



E 



\TAn,K)\'^ (4) 



X:sup 



across finer scales 



where sum takes supremum over so-called maxima lines, 
i.e., lines which link the local maxima at the present scale 
with the maxima at the previous scales. 

— In Multifractal Detrended Fluctuation Analy- 
sis (MDFA) approach a multiresolution coefficient is de- 
fined as, [8|]- 



SpK{n,K)=i^[Xii)-P^{i)f 



and then the partition function 6'(n, q) is given as 
F{n,q)=J2SpK{n,Ky 

K 



(5) 



(6) 



Let us notice that the supremum over the maxima 
lines, calculated in WTMM approach (|4]), provides the 
inter scale links while MDFA method works on each scale 
independently. 

Dependence on n of any partition function P leads to 
so-called structure function r^q) as follows: 



P{n,q) 



and then by applying a Legendre transform to a curve 
(9) oiie obtains the multifractal spectrum {h, D{h)): 



h- 



dTjq) 
dq 



D{h) =qh- T{q) 



(8) 



which links the singularity exponent value h with its 
probability to appear D{h) in a signal. 

In case of a monofractal signal one expects a linear 
structure function r(g) and then a corresponding spec- 
trum should consists of one point (/i, D{h)). But the nu- 
merics overestimate the real Holder exponent h. There- 
fore the numerical results are sometimes called the coarse 
Holder exponents [4]. Additionally, the Legendre trans- 
form (|8]) significantly magnifies any small departures 
from linearity. Finally, what we obtain from numerical 
estimates are spectra which are far from being a point- 
like pLO, 20]. However the maximum of a spectrum curve 
is directly related to the theoretical spectrum point and 
the maximum point is the highest occupied point of the 
whole spectrum [20]. A multifractal spectrum is distinct 
from a monofractal plot by existence of many highly oc- 
cupied points. This feature is called the linearization 
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FIG. 4: Structure functions r(^) obtained by MDFA for 
a RR signal. (A) T{q) for nsr_gda03_wake. (B) t(^) 
for nsr_gda03_wake shuffled at random. (C) T{q) for 
nsr_gda03_wake integrated. (D) T{q) for nsr_gda03_wake shuf- 
fled at random and integrated. All linear fits provides r^- 
errors greater than 0.98 



effect and is suggested to be related with the nature of 
the process Q. 

In general when studying empirical data what we know 
is that negative moments, i.e., when (7 < 0, are respon- 
sible for enhancing small fluctuations while positive mo- 
ments, <7 > 0, amplify large fluctuations. The right wing 
of a spectrum corresponds to negative moments and the 
left wing reflects properties collected for positive mo- 
ments. 

In FigHlA we present the structure function found due 
to MDFA method to a series shown in Fig|2l The struc- 
ture function T{q) of the pure signal is concave and any 
linear fit is loaded by a large -error. However, one can 
easily see that for > 4 line fits are justified (the lin- 
earization effect) what will result in emergence of two 
dense occupied points on (h^D(h)) plane: at = 0.290 
and at = 0.005. The signal is multifractal. If the data 
is shuffled at random, then the corresponding partition 
function becomes a flat line, see Fig. [4^. The same result 
is obtained for a white noise — a canonical example of a 
monofractal signal. Its multifractal spectrum consists of 
a point (0, 1). 

In the original detrended fluctuation analysis one deals 
with so-called profiles, see The profile can be 

thought as a signal obtained after integrating series. 
Dealing with series of integrated RR intervals one intro- 
duces a kind of smoothing since only growth in a data set 
is possible. Moreover, integrated series allows to study 
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strength of independency in a series. In Fig. [4lG the 
structure function obtained for the integrated series of 
nsr_gda03_wake is shown. Again the hnear fits are pro- 
vided for 1^1 > 4. The hnear approximations in these 
intervals give: 1.227 and 1.019 what locates the multi- 
fractal spectrum between these two points. Hence, the 
spectrum is moved right by about of 1 from the spectrum 
obtained from the pure series on the (h^D{h)) plane. 
However, if the data is independent then, see the struc- 
ture function of integrated shuffled data in Fig. HP, the 
linear fit provides the slope of 0.49 what means that the 
spectrum is moved right but only by 1/2. 



IV. HEART INTERBEAT INTERVAL DATA 

The series called nsr_gda were collected from 39 
healthy individuals (4 women, 35 men, the average age 
is 54 zb 7) without past history of cardiovascular disease, 
with both echocardiogram and electrocardiogram in nor- 
mal range. The left ventricle ejection fraction was nor- 
mal (mean LVEF = 68.0 ± 4.7%). For each person the 
24-hour EGG Holter monitoring was performed. The sig- 
nal was digitized using Delmar Avionics recorder (Digi- 
torder) and then analyzed and annotated using Delmar 
Accuplus 363 system (fully interactive method) by an ex- 
perienced physician. The minimum number of qualified 
normal sinus beats required for the signal to enter the 
study is 85%. Moreover, the standard filters [l^ have 
been applied to smooth out the data. 

From each signal the two 5-hour continuous subsets 
were extracted: diurnal and nocturnal. By extracting 
such subseries we wanted eliminate at least some sources 
of nonstationarities. The subsets were extracted manu- 
ally, with regards to two easily recognizable stages: daily 
activity and sleep state. In most cases these periods are 
consistent with hours 15:00-20:00 for the daily activity 
and 0:00-5:00 for the sleep state. In the following we 
will denote the subseries as nsr_gda_wake, nsr_gda_sleep, 
respectively. All RR-series and their subseries are acces- 
sible on request [28]. 

Since the WTMM method requires data sets longer 
than 2^^, then some of the sleep subseries have to be ex- 
cluded from WTMM calculations. So that nsr_gda_sleep 
group consists of 34 data sets in case of WTMM calcula- 
tions. 



V. SCALING PROPERTIES OF PARTITION 
FUNCTIONS 

Both partition functions Z{n^ q) and F(n, q) are calcu- 
lated by tools accessible from PhysioNet [23] • Namely, we 
work with two Physionet packages: mult ifr act al.c (pre- 
pared by Y. Ashkenazy) and dfa.c (prepared by J.Mietus, 
G-K Peng, and G. Moody ) adapted by us and accessible 
from [2a] 





FIG. 5: Group average of partition functions calculated by 
MDFA and WTMM methods for each element of the group 
and for different q = —5,— 4,. ..,4, 5. In green, statistical 
errors of estimates (standard deviation of a group \^/N, where 
N number of elements in the group) are given. 



For each individual series, regarded as either directly as 
a sequence of RR intervals, or as a sequence of integrated 
RR intervals, the partition function was found for all q in 
[—5, 5] with a step 0.1, for scales from n = 10 to almost 
1000 heart beats. In all calculations we applied the third 
derivative of a Gaussian as a mother wavelet in WTMM 
method, (thus up to quadratic polynomial dependencies 
in t are eliminated by WTMM) and a quadratic poly- 
nomial as a detrending fit (MDFA). Such setting of pa- 
rameters has been found by us as less dependent on the 
length of test series — localization of spectra is stable 
(WTMM) and width of spectra is stable (MDFA). 

The partition functions Z(n, g) and F(n, g), shown in 
Figs mini are group averages obtained from the partition 
functions found for each individual series by WTMM and 
MDFA methods, respectively. We present Z{n^ q) and 
F(n, q) for ^'-representatives {—5, —4, . . . , 4, 5} together 
with their statistical errors. We found errors satisfactory 
low (with significance level p < 0.05 ) to discriminate val- 
ues of particular gs. Therefore we could turn to searching 
for power-law scalings of these functions. To reliably in- 
fer power-law scaling, a straight line in log-log plot has 
to be established. In Figs [5l [6] the log- log scale is used to 
enhance visually linear relations between log [scale] and 
log [partition function]. 

It appears that the partition functions F(n, q) change 
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FIG. 6: Group average of partition functions calculated by 
MDFA and WTMM methods for each integrated signal of the 
group for different 5,— 4,...,4,5. In green, statistical 

errors of estimates (standard deviation of a group \\^, where 
N number of elements in the group) are given. 
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their dependence on a scale. Especially, the change is 
noticeable in case of pure series. Since numerical cal- 
culations are sensitive to details of curves, we evaluate 
subsequently local slopes of all functions: logF(n,g') vs 
logn as well as log Z{n^ q) vs logn to observe changes in 
scalings. In particular, we calculate linear fits for sub- 
sequent 20 points of each partition function and move 
with the starting point from 20th heart beat to the end 
of series. This way we obtain the variety of possible ap- 
proximations and then we can ask how one can estimate 
fine scale or/and large scale limits properties. 

The results - local structure functions and resulting 
local spectra plots are shown in Figs [71 [HI The first struc- 
ture function TfirstiQ)^ labeled first (red plots) starts at 
the 20th heart beat. The plots labeled last in the figures 
correspond to riast{Q) obtained in the last possible but 
smaller than 1000th heart beat (green plots). One has 
to be aware that 20 points at fine scales correspond to 
about several heart beats while 20 points at large scales 
cover few hundred heart beats since a scale box grows 
geometrically with a ratio 2^-^^. The spectra are plotted 
below the corresponding structure functions. One can see 
that the 'body' of each spectrum occupies a wide part of 
a plane (h^D(h)). The strange structures appearing in 
some of spectra are consequences of Legendre transfor- 
mation ([8|) when treated by numerical methods. 

Both methods MDFA and WTMM provide similar re- 



FIG. 7: The structure functions and spectra obtained when 
the linear fits are found for subsequent 20 points of the corre- 
sponding partition functions. Group averages for integrated 
nsr_gda_wake obtained by both methods for pure series (upper 
panel) and integrated series (bottom panel). 

suits in case of pure series. The fine scales are signifi- 
cantly separated from large scales: the first-plots are cen- 
tered at = 0.6 while last-plots are concentrated around 
h = 0.2. The maxima of corresponding spectra are lo- 
cated similarly on {h^D{h)) -plane and widths of corre- 
sponding spectra are in similar relations: the fine scale 
widths are significantly greater than the large scale spec- 
tra. These results might suggest presence of plenty in- 
dependent stochastic noises actively influencing dynamic 
properties of heart rate in short time scales. 

The integrated series are expected to provide spectra 
moved right by 1 from the pure series. It appears that 
this simple relation is about to be satisfied only in case 
of the spectra of wake series treated by MDFA method. 
In case of the sleep data we see that the spectrum cor- 
responding to large scales is located in h < 1 area. This 
observation suggests that nocturnal series are character- 
ized by stronger independence between subsequent heart 
beats than diurnal series. In case of WTMM method 
both limits: fine and large scale provide spectra centered 
similarly at about h = 1.2. Therefore one can conjecture 
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FIG. 8: The structure functions and spectra obtained when 
the Hnear fits are found for subsequent 20 points of the corre- 
sponding partition functions. Group averages for integrated 
nsr_gda_wake obtained by both methods for pure series (upper 
panel) and integrated series (bottom panel). 



that WTMM spectra of integrated series represent some 
average mult ifr act al characteristics. 



VI. SPECTRA 

In the following we search for the best linear fits of 
partition functions in the particular time intervals. These 
intervals, transfered into number of beats, are directly 
related to the physiology of the heart rate control and 
the classical spectral analysis as it was described in Sec 
II. The group structure functions and the group spectrum 
plots are presented as follows, see FiglQl 

• 10-25 heart beats that cover times from 8 to 20 
seconds, what allows us to estimate early compen- 
satory mechanism represented by LF. 

• 25-300 heart beats that last from 20 seconds to 
4 minutes, what examines the activity of rennin- 
angiotensin-aldosterone system and other neuro- 



hormones, which additionally activates the sympa- 
thetic nervous system and corresponds to VLF. 

• over 300 heart beats, i.e., times over 4 minutes allow 
to study frequencies smaller than 0.003Hz, hence 
we are within ULF component. 

We can not enter into the high frequency HF range 
because of numerical instabilities occurring when q < 0. 

It appears that in short times, less than 20 seconds, 
each method locates the multifractal spectra in different 
intervals of h. Both pure spectra: wake and sleep are in 
h e (0.5, 1.0) in case of MDFA while maxima of WTMM 
spectra are about h = 0.2 though there is a noticeable 
center of WTMM wake spectrum points at h = 0.6. The 
spectra obtained from integrated series sustain this obser- 
vations — WTMM plots are significantly separated from 
plots obtained by means of MDFA. It denotes that rough- 
ness of signals measured in scales smaller than 25 beats 
(a half of a box from Fig|3]) is considerable differently 
estimated by both methods. Yet both methods indicates 
on presence of strong persistence in series, because the 
spectra contain singularity exponents with h > 0.5. The 
persistence is usually related to burstiness or clustering 
in data — fluctuations larger or smaller than a given level 
happen subsequently, i.e., much often than in the case of 
independent observations. 

Heart rate in LF range is strongly influenced by os- 
cillations in blood pressure. Basically the baroreflex is 
the modulator of RR interval in a healthy subject. How- 
ever there are receptors which control the other aspects 
of proper blood supply, and which activity can be evalu- 
ated by the level of randomness in LF range [l^. Both 
multifractal methods produce wide spectra what can be 
attributed to a large variety of controlling elements under 
which a cardiovascular system maintains blood supply 
during daily human activity. Moreover, the randomness 
point h = 0.5 belongs to both spectra and it could pro- 
vide an indicator of the level of the baroreflex sensitivity. 
MDFA plots of integrated series are moved toward a mul- 
tifractal range which corresponds to a series of integrated 
random walk. Thus RR intervals are seen as rather inde- 
pendent of each other. However WTMM provides plots 
with maxima highly occupied at about = 0.85 for sleep 
series and around = 0.95 in case of wake series. 

DFA picture of VLF range, Fig|9]plots of middle panel, 
again provides a huge variety of scaling exponents in case 
of daily activity h G (0, 0.6) with two centers where spec- 
trum points are densely located: at maximum h = 0.3 
and at h = 0.6. WTMM spectrum in VLF range gives 
the wake spectrum narrower and mainly concentrated at 
its left wing, i.e., for h < 0.2. Thus comparing the two 
methods one can say that WTMM enhances influence 
of large oscillations which are contained in r{q) when 
g > 0, while MDFA amplifies role of small fluctuations 
expressed by r{q) with q < 0. Both wake spectra are 
wider than the spectra of sleep series and the sleep spec- 
tra are located inside the corresponding wake spectra. 
Thus one can conjecture loosing variability in heart rate 
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FIG. 9: The structure functions and spectra obtained when 
the Unear fits are for heart beat intervals physiologically mo- 
tivated. The intervals are described in plot titles. Every third 
points are plotted to emphasize centers of spectral point con- 
centration. 



(measured by a pool of scaling exponents) during the 
night rest in VLF range. 

The integrated wake and sleep spectra corresponding 
to VLF component exhibit noticeable separation from 
each other: the sleep spectra are moved to left, to lower 
Holder exponents, namely, h G (0.6, 1). Both sleep curves 



are characterized by densely occupied left wings. Hence 
the roughness of night series is significantly wilder. The 
main parts of the wake spectra occupy the similar range 
of h, i.e., h e (1,1.2). This means that WTMM wake 
spectrum is moved right by 1 when compared it to the 
pure series result but MDFA one not. However there is 
a spectrum center at /i = 1.5. 

The third frequency component UHF corresponds to 
the time range where long range dependences are usually 
searched. Therefore it is not a surprise that our results 
recalls results obtained by others, see [13, E, 0, US] • 
case of pure series all spectra are characterized by similar 
set of singularity exponents h concentrated at h = 0.2 
and the widths of WTMM spectra are wider than the 
corresponding MDFA ones. Moreover, MDFA results re- 
ceived for integrated series reveal global properties known 
thanks to DFA studies [25[ : the healthy heart long range 
correlations are characterized hy a^FA = ^ ^oi di daily 
activity and by aoFA = 0.8 during a nigh rest. WTMM 
spectra for integrated series differentiates day and night 
series by spectral width: the sleep spectrum is signifi- 
cantly wider than the wake one. However, the left wings 
of these spectra - consisting of many points, are located 
at the same h as it is found by DFA. 



VII. CONCLUSIONS 

Many naturally occurring empirical time series violate 
the assumption of independence of subsequent elements. 
Instead the data often point at the presence of strong 
temporal correlations. Relating these correlations with 
dynamic properties of a system under study is challeng- 
ing. 

Heart rate variability is a measure of fundamental as- 
pects of a human physiology. It is known that frequency 
domain analysis of heart rate provides predictors of mor- 
tality in chronic heart failure — LF is a predictor for 
sudden cardiac death while VLF and ULV are predictors 
for all-causes mortality, see jl], 05 E] • 

In the following, basing on the descriptive rather as- 
pect of multiscaling, we have analysed scaling properties 
of average multifractal partition functions in these phys- 
iologically grounded interbeat intervals: LF, VLF and 
ULF, following normal RR intervals of 39 healthy sub- 
jects. The data was grouped according to the circadian 
rhythm (wake and sleep series), the methods of multi- 
fractal analysis (WTMM and MDFA approach), and RR 
intervals were compared to integrated RR intervals. By 
applying different numerical methods we wanted to de- 
termine some objective properties hidden in series con- 
sidered which discriminate daily heart activity from its 
nocturnal rest. In general, we observed loss of heart rate 
variability during night in VLF range and its increase in 
time corresponding to ULF component. 

However, we have also found that details of multi- 
fractal properties such as spectrum width and location 
are different for different numerical methods. A parti- 
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tion function obtained by WTMM remembers about the 
strength of oscihations in a signal read in finer scales 
while MDFA provides at each scale an independent par- 
tition function. Therefore the relations between scales 
are more consistent with each other in case of WTMM 
spectra seen in different scales. It is specially noticeable 
in case of integrated series. The study of integrated RR- 
intervals has put a significant impact in evaluation of the 
participation of independent stochastic sources. We hope 



that all these differences would lead to additional tools 
allowing to quantify and qualify neural and humoral con- 
trol of heart rate. 
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